<html>
<head>
<script src="~/three.js/build/three.min.js"></script>
<script src="~/three.js/examples/js/controls/TrackballControls.js"></script>

<script type="text/python">
pythonjs.configure(javascript=True)

ren = None
scn = None
cam = None

def avg( arr ):
	n = arr.length
	x = []; y = []; z = []
	for v in arr:
		x.push( v.x )
		y.push( v.y )
		z.push( v.z )
	return new THREE.Vector3( sum(x)/n, sum(y)/n, sum(z)/n )

def avg2( a, b ):
	x = a.x + b.x
	y = a.y + b.y
	z = a.z + b.z
	return new THREE.Vector3(x/2, y/2, z/2)

def avg4( a, b, c, d ):
	x = a.x + b.x + c.x + d.x
	y = a.y + b.y + c.y + d.y
	z = a.z + b.z + c.z + d.z
	return new THREE.Vector3(x/4, y/4, z/4)

class CatmullClark:
	'''
	http://www.rorydriscoll.com/2008/08/01/catmull-clark-subdivision-the-basics/

	'''

	def id(self, a,b):
		v1 = Math.min(a,b)
		v2 = Math.max(a,b)
		return v1+':'+v2

	def __init__(self, geo, parent=None):
		self.input_geom = geo
		geo.computeCentroids()

		# Q is the average of the surrounding face points
		q = [ [] for i in range(geo.vertices.length) ]

		# R is the average of the surrounding edge midpoints
		r = [ [] for i in range(geo.vertices.length) ]

		self.edge_faces = edge_faces = {}
		for face in geo.faces:
			id1 = self.id(face.a, face.b)
			id2 = self.id(face.b, face.c)
			id3 = self.id(face.c, face.a)
			if id1 not in edge_faces:
				edge_faces[id1] = []
			if id2 not in edge_faces:
				edge_faces[id2] = []
			if id3 not in edge_faces:
				edge_faces[id3] = []

			edge_faces[id1].append( face )
			edge_faces[id2].append( face )
			edge_faces[id3].append( face )

			mid = avg2( geo.vertices[face.a], geo.vertices[face.b] )
			r[ face.a ].push( mid )
			r[ face.b ].push( mid )

			mid = avg2( geo.vertices[face.b], geo.vertices[face.c] )
			r[ face.b ].push( mid )
			r[ face.c ].push( mid )

			mid = avg2( geo.vertices[face.c], geo.vertices[face.a] )
			r[ face.c ].push( mid )
			r[ face.a ].push( mid )

			q[ face.a ].push( face.centroid )
			q[ face.b ].push( face.centroid )
			q[ face.c ].push( face.centroid )


		# 0.5  from the surrounding edge midpoints `R`, 
		R = []
		for a in r:
			R.push( avg(a).multiplyScalar(0.5) )

		# 0.25 weight from the average of surrounding face-points `Q`, 
		Q = []
		for u in q:
			Q.push( avg(u).multiplyScalar(0.25) )

		fgeo = new THREE.Geometry()
		egeo = new THREE.Geometry()
		rgeo = new THREE.Geometry()

		sgeo = new THREE.Geometry()

		edgeverts = []
		i = 0
		for face in geo.faces:
			id = self.id(face.a, face.b); edge = edge_faces[id]
			e1 = avg4( geo.vertices[face.a], geo.vertices[face.b], edge[0].centroid, edge[1].centroid )

			id = self.id(face.b,face.c); edge = edge_faces[id]
			e2 = avg4( geo.vertices[face.b], geo.vertices[face.c], edge[0].centroid, edge[1].centroid )

			id = self.id(face.c, face.a); edge = edge_faces[id]
			e3 = avg4( geo.vertices[face.c], geo.vertices[face.a], edge[0].centroid, edge[1].centroid )

			egeo.vertices.append(e1)
			egeo.vertices.append(e2)
			egeo.vertices.append(e3)

			fgeo.vertices.append( face.centroid )

			sgeo.vertices.append(e1)
			sgeo.vertices.append(e2)
			sgeo.vertices.append(e3)
			sgeo.faces.push( new THREE.Face3(i,i+1,i+2) )


			# 0.25 from the original control-point.

			vpoint = R[face.a].clone().add( Q[face.a] ).clone().add( geo.vertices[face.a].clone().multiplyScalar(0.25) )
			sgeo.vertices.push( vpoint )
			rgeo.vertices.push( vpoint )

			vpoint = R[face.b].clone().add( Q[face.b] ).clone().add( geo.vertices[face.b].clone().multiplyScalar(0.25) )
			sgeo.vertices.push( vpoint )
			rgeo.vertices.push( vpoint )

			vpoint = R[face.c].clone().add( Q[face.c] ).clone().add( geo.vertices[face.c].clone().multiplyScalar(0.25) )
			sgeo.vertices.push( vpoint )
			rgeo.vertices.push( vpoint )

			sgeo.faces.push( new THREE.Face3(i+3,i,i+2) )
			sgeo.faces.push( new THREE.Face3(i+4,i+1,i) )
			sgeo.faces.push( new THREE.Face3(i+5,i+2,i+1) )


			i += 6



		pmat = new THREE.ParticleSystemMaterial( size=1, color=0x0000ff )
		p = new THREE.ParticleSystem( egeo, pmat )
		parent.add( p )

		pmat = new THREE.ParticleSystemMaterial( size=0.5, color=0xffff00 )
		p = new THREE.ParticleSystem( rgeo, pmat )
		parent.add( p )

		pmat = new THREE.ParticleSystemMaterial( size=0.4, color=0xff0000 )
		p = new THREE.ParticleSystem( fgeo, pmat )
		parent.add( p )

		sgeo.mergeVertices()
		#sgeo.computeCentroids();
		sgeo.computeFaceNormals();
		sgeo.computeVertexNormals()
		mat = new THREE.MeshLambertMaterial( color=0x00ff00, wireframe=false, wireframeLinewidth=3 )
		self.mesh = new THREE.Mesh( sgeo, mat )
		parent.add( self.mesh )
		self.output_geom = sgeo
		self.mesh = mesh

def main():
	global ren, scn, cam, mesh, controls

	div = document.createElement( 'div' )
	document.body.appendChild(div)

	width = window.innerWidth-20; height = window.innerHeight-20
	scn = new( THREE.Scene() )
	cam = new( THREE.PerspectiveCamera( 45, width/height, 0.01, 10000) )
	cam.position.z = 60
	cam.position.x = 5
	controls = new THREE.TrackballControls( cam )

	ren = new( THREE.WebGLRenderer() )
	ren.setSize( width, height )

	div.appendChild( ren.domElement )

	light = new( THREE.PointLight(0xffffff, 2, 500) )
	light.position.set( 0, 100, 90 )
	scn.add( light )

	geo = new THREE.BoxGeometry( 20,20,20 )
	mat = new THREE.MeshBasicMaterial( wireframe=True )
	mesh = new THREE.Mesh( geo, mat )
	scn.add( mesh )

	ss1 = CatmullClark( geo, parent=mesh )
	ss2 = CatmullClark( ss1.output_geom, parent=scn )
	ss2.mesh.position.x += 30

	animate()

def animate():
	requestAnimationFrame( animate )
	controls.update()
	ren.render( scn, cam )

</script>

</head>
<body onload="main()">
</body>
</html>